% Ancillary script for "The Climate Risk Premium" (Lemoine, JAERE, 2020).

% Calculates average trajectories of variables as well as their quantiles.
% Also fits lognormal distribution to integral of temperature over first 50
% years, used in calibrating damages.

% Begin from loaded workspace

% For Figure 2, that workspace should be created with consumption volatility only, with option_otherloop=montecarlo_damage, or with option_otherloop=montecarlo_cs

% For damages fitting, run with option_otherloop=montecarlo_cs and option_nodamages='on'


%% Distribution of key variables over time

% Take expectations over damages, so obtain matrices with dimensions simulation x time x climsens_node
temp_co2 = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 length(weights_damage)]),Paths.M1(:,:,:,:)+Paths.M2(:,:,:,:)),4 );
temp_temp = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 length(weights_damage)]),Paths.temp(:,:,:,:)),4 );
temp_cons = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 length(weights_damage)]),Paths.cons(:,:,:,:)),4 );

% Take expectations over climate sensitivity, so obtain matrices with dimensions simulation x time
temp_co2 = sum( bsxfun(@times,reshape(weights_climsens,[1 1 length(weights_climsens)]),temp_co2(:,:,:)),3 );
temp_temp = sum( bsxfun(@times,reshape(weights_climsens,[1 1 length(weights_climsens)]),temp_temp(:,:,:)),3 );
temp_cons = sum( bsxfun(@times,reshape(weights_climsens,[1 1 length(weights_climsens)]),temp_cons(:,:,:)),3 );

quantile_low = 0.25;
quantile_high = 0.75;

% Calculate statistics at each timestep
mean_co2(:,1) = mean(temp_co2,1)'/Params.gtc_per_ppm; % ppm
quantile_co2(:,1) = quantile(temp_co2,quantile_low,1)'/Params.gtc_per_ppm; % ppm
quantile_co2(:,2) = quantile(temp_co2,quantile_high,1)'/Params.gtc_per_ppm; % ppm
quantile_co2(:,3) = min(temp_co2,[],1)'/Params.gtc_per_ppm; % ppm
quantile_co2(:,4) = max(temp_co2,[],1)'/Params.gtc_per_ppm; % ppm
mean_temp(:,1) = mean(temp_temp,1)'; % deg C
quantile_temp(:,1) = quantile(temp_temp,quantile_low,1)'; % deg C
quantile_temp(:,2) = quantile(temp_temp,quantile_high,1)'; % deg C
quantile_temp(:,3) = min(temp_temp,[],1)'; % deg C
quantile_temp(:,4) = max(temp_temp,[],1)'; % deg C
mean_cons(:,1) = mean(temp_cons,1)'/10^12; % in trillion dollars
quantile_cons(:,1) = quantile(temp_cons,quantile_low,1)'/10^12; % in trillion dollars
quantile_cons(:,2) = quantile(temp_cons,quantile_high,1)'/10^12; % in trillion dollars
quantile_cons(:,3) = min(temp_cons,[],1)'/10^12; % in trillion dollars
quantile_cons(:,4) = max(temp_cons,[],1)'/10^12; % in trillion dollars
mean_cons_pc(:,1) = mean_cons(:,1)*10^(12)./transpose(Paths.pop)/10^3; % in thousand dollars p.c.
quantile_cons_pc(:,:) = bsxfun(@rdivide, quantile_cons(:,:)*10^(12)/10^3, transpose(Paths.pop)); % in thousand dollars p.c.

clear temp_co2 temp_cons;


%% Now fit a lognormal distribution to the integral of temperature over the first 50 years

integral_temp = Params.timestep*trapz(temp_temp(:,[1:50/Params.timestep]),2); % simulations x 1
log_integral_temp = log(integral_temp);
[fitted_normal.mu,fitted_normal.sigma,fitted_normal.mu_ci,fitted_normal.sigma_ci] = normfit(log_integral_temp);

alpha_upperbound = -log(0.5)/mean(integral_temp);


% % plot fit
% samples_log_integral_temp = normrnd(fitted_normal.mu,fitted_normal.sigma,Params.montecarlo_draws,1);
% histogram(log_integral_temp)
% hold on;
% histogram(samples_log_integral_temp)
% legend('actual','fitted');
% xlabel('log integral temp');
% ylabel('count')

%clear temp_temp;


%% Fit level damages

% Sample from Pindyck (2019) lognormal distribution for his damage
% parameter phi, truncating at same point used for main text
for index_sample=1:Params.montecarlo_draws
    [phi_sample(index_sample,1), seed] = truncated_normal_b_sample( -2.446, 1.476, log(-log(0.5)), 123); % need use original scale but adjusted location parameter
end
phi_sample = exp(phi_sample); %  adjust for having been lognormal
bT2_sample = 1-exp(-phi_sample);

T2_sample = temp_temp(:,50/Params.timestep).^2;

[b_fitted_lognormal.mu,b_fitted_lognormal.sigma,b_fitted_lognormal.mu_ci,b_fitted_lognormal.sigma_ci] = normfit(log(bT2_sample./T2_sample));

% find value of b corresponding to upper bound of 50% consumption loss at mean
% of T^2
b_upperbound = 0.5/mean(T2_sample);

% find ranking of Nordhaus parameter
[~,idx]=ismember(0.00236,sort([0.00236; bT2_sample./T2_sample],'ascend'));
nordhaus_pctile = idx/(1+length(bT2_sample));

% % plot fit
% samples_log_b = normrnd(b_fitted_lognormal.mu,b_fitted_lognormal.sigma,Params.montecarlo_draws,1);
% histogram(log(bT2_sample./T2_sample))
% hold on;
% histogram(samples_log_b)
% legend('actual','fitted');
% xlabel('ln(b)');
% ylabel('count')

% % plot result
% samples_log_b = normrnd(b_fitted_lognormal.mu,b_fitted_lognormal.sigma,Params.montecarlo_draws,1);
% histogram(exp(samples_log_b))
% xlabel('b');
% ylabel('count')
% xlim([0 0.1])


% % plot original pindyck distribution
% histogram(1-exp(-phi_sample))
% xlabel('Fractional consumption loss in 2050');
% ylabel('count')


%% Now calculate the expected trajectories and standard deviations at each damage quadrature node

% Take expectations over climate sensitivity, so obtain matrices with
% dimensions simulation x time x damage_node
temp_co2 = sum( bsxfun(@times,reshape(weights_climsens,[1 1 length(weights_climsens) 1]),Paths.M1(:,:,:,:)+Paths.M2(:,:,:,:)),3 );
temp_temp = sum( bsxfun(@times,reshape(weights_climsens,[1 1 length(weights_climsens) 1]),Paths.temp(:,:,:,:)),3 );
temp_cons = sum( bsxfun(@times,reshape(weights_climsens,[1 1 length(weights_climsens) 1]),Paths.cons(:,:,:,:)),3 );

% Calculate statistics at each timestep, for each damage node
mean_co2_damnode = squeeze(mean(temp_co2,1))/Params.gtc_per_ppm; % ppm
std_co2_damnode = squeeze(std(temp_co2,0,1))/Params.gtc_per_ppm; % ppm
mean_temp_damnode = squeeze(mean(temp_temp,1)); % deg C
std_temp_damnode = squeeze(std(temp_temp,0,1)); % deg C
mean_cons_damnode = squeeze(mean(temp_cons,1))/10^12; % in trillion dollars
std_cons_damnode = squeeze(std(temp_cons,0,1))/10^12; % in trillion dollars
mean_cons_pc_damnode = bsxfun(@rdivide,mean_cons_damnode,transpose(Paths.pop))*10^(12)/10^3; % in thousand dollars p.c.
std_cons_pc_damnode = bsxfun(@rdivide,std_cons_damnode,transpose(Paths.pop))*10^(12)/10^3; % in thousand dollars p.c.
    
clear temp_co2 temp_temp temp_cons;


%% Now calculate the expected trajectories and standard deviations at each climate sensitivity quadrature node

% Take expectations over damages, so obtain matrices with
% dimensions simulation x time x cs_node
temp_co2 = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 length(weights_damage)]),Paths.M1(:,:,:,:)+Paths.M2(:,:,:,:)),4 );
temp_temp = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 length(weights_damage)]),Paths.temp(:,:,:,:)),4 );
temp_cons = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 length(weights_damage)]),Paths.cons(:,:,:,:)),4 );

% Calculate statistics at each timestep, for each cs node
mean_co2_csnode = squeeze(mean(temp_co2,1))/Params.gtc_per_ppm; % ppm
std_co2_csnode = squeeze(std(temp_co2,0,1))/Params.gtc_per_ppm; % ppm
mean_temp_csnode = squeeze(mean(temp_temp,1)); % deg C
std_temp_csnode = squeeze(std(temp_temp,0,1)); % deg C
mean_cons_csnode = squeeze(mean(temp_cons,1))/10^12; % in trillion dollars
std_cons_csnode = squeeze(std(temp_cons,0,1))/10^12; % in trillion dollars
mean_cons_pc_csnode = bsxfun(@rdivide,mean_cons_csnode,transpose(Paths.pop(1,[1:size(Paths.cons,2)])))*10^(12)/10^3; % in thousand dollars p.c.
std_cons_pc_csnode = bsxfun(@rdivide,std_cons_csnode,transpose(Paths.pop(1,[1:size(Paths.cons,2)])))*10^(12)/10^3; % in thousand dollars p.c.
    
clear temp_co2 temp_temp temp_cons;


